Liquid-liquid coexistence in the phase diagram of a fluid confined 

in fractal porous materials. 
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Multicanonical ensemble sampling simulations have been performed to calculate the phase diagram 
of a Lennard- Jones fluid embedded in a fractal random matrix generated through diffusion limited 
cluster aggregation. The study of the system at increasing size and constant porosity shows that 
the results are independent from the matrix realization but not from the size effects. A gas-liquid 
transition shifted with respect to bulk is found. On growing the size of the system on the high 
density side of the gas-liquid coexistence curve it appears a second coexistence region between two 
liquid phases. These two phases are characterized by a different behaviour of the local density inside 
the interconnected porous structure at the same temperature and chemical potential. 

PACS numbers: 61.20.Ja, 61.20.-p, 64.70.Fx 



I. INTRODUCTION 

The effect of confinement on the phase diagram of flu- 
ids is a longstanding problem with a wide interest from 
the theoretical point of view and a wide area of possible 
applications such as catalysis, adsorption and filtration. 
A large phenomenology on this subject can be foundi. 
Real porous materials show generally a complex struc- 
ture made of an interconnected network of pores. Vycor 
glass and silica gels are relevant examples. It is impor- 
tant to distinguish between two classes of systems. In 
hosting media with very high porosity like silica aerogel 
(90% — 98%) it has experimentally been observed that 
the gas-liquid transition is preserved although the con- 
finement causes a reduction of the critical temperature 
and a shrinkage of the gas-liquid coexistence curveS^S'. In 
mesoporous materials with a porosity in the range be- 
tween 30% and 60%, as for instance Vycor glass, there 
are not direct observations of equilibrium phase transi- 
tions^. This last behaviour is predicted by theoretical 
work performed in the framework of lattice gas mod- 
gj r)i,4,5. 6.7.8.9^ rpj^jg type of approach is connected to the 
more general issue of the effects of quenched disorder on 
critical phenomenaiSiiiiiSiiii^. In dilute silica aerogels 
with high porosity, however, phase transitions between 
equilibrium phases cannot be theoretically excludedi^ii^ 
although more recent experimentsiii^ and mean field 
calculations-''^ questioned the early experimental findings 
of a gas-liquid transition in ^He in aerogel. 

For dilute silica aerogels off lattice liquid state mod- 
els have been introducedii'^° to take into account more 
in detail the microstructure of the confining systems. In 
these off lattice models the disordered porous material 
is modeled as a system of spheres frozen in a prede- 
fined structure. The quenched matrix approximately ac- 
counts for the geometric constraints of the interconnected 
pore structure of the solid. From both computer simu- 
lationsSfliSi and integral equation methodsi^ it has been 
found that the phase diagram of these models can show, 



besides the gas-liquid coexistence (GLC), the presence of 
a second transition. This transition can be a gas-gas or a 
liquid-liquid coexistence depending on the liquid-matrix 
type of interaction. Nonetheless the appearance of this 
coexistence and its details depend on the way in which 
the confining random matrix is realized^^. Moreover rel- 
evant wetting or drying effects could be present^S. 

Therefore while the coexistence of different fluid forms 
in a simple monatomic fluid confined in silica aerogels is 
of uttermost interest, its existence is still a controversial 
result. The fact that the phase diagram has shown to 
depend so much on the details of the quenched random 
matrix has precluded so far more detailed microscopic 
analysis on these confined fluids and a systematic study 
of size effects. 

From these considerations and a more recent theo- 
retical analysis based on integral equation methods^ it 
emerges the importance of realizing by computer simu- 
lation an host system which accounts microscopically for 
the high porosity and the fractal structure of the silica 
aerogel. Its structure is due to the process of formation 
of the gel through the random aggregation of the sil- 
ica particles and it is measured in small angle scattering 
experiments where a fractal dimension oi D = 1.8 has 
been estimated22i2^*2^. These microscopic features are 
not present in any of the previous lattice or off lattice 
simulated systems, mentioned above. 

In order to build up a confining medium with the struc- 
ture of silica aerogels, we have implemented a numerical 
procedure^^ based on the diffusion hmited cluster-cluster 
aggregation (DLCA)SiS^i2^. It has been shown2i2i2^ 
that the DLCA algorithm is able to reproduce the struc- 
ture factor and the fractal dimension of silica aerogels. 
In a previous study we calculated with the Multicanoni- 
cal Ensemble Samphng (MES)2LSS*S the phase diagram 
of a Lennard-Jones fluid confined in the matrix gener- 
ated with the DLCA and we found a GLC curve2& and 
no evidence of a second transition except for a possible 
signature in a shouldering on the liquid side of the coex- 
istence curve of a liquid-liquid coexistence (LLC). As it 
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is well known for aerogels^* systems with the same poros- 
ity show the same average size of the clusters, the same 
connectivity, the same behaviour and range of the fractal 
scale. We expect therefore that performing simulations 
with confining media at constant porosity and increasing 
size that the behaviour of the liquid will be independent 
from the details of the realization of the matrix apart for 
the size effects themselves since the liquid is embedded 
in self similar structures. 

In this paper we report the results obtained by increas- 
ing the size of the confining medium at constant porosity 
with the aim of inquiring about the existence of a LLC 
in the Lennard- Jones simple fluid possibly hidden in our 
past study by the small size of the simulated system. 



II. SIMULATION DETAILS 

To generate the off lattice matrix in a simulation box of 
side L we introduce a number Ns of spheres of diameter 
(Ts randomly placed and apply periodic boundary condi- 
tions. After the application of the DLCA algorithm the 
final configuration of the matrix consists of a percolated 
cluster of Ng spheres. For our former system, which we 
will call in the following system we used L = IStJs 
and Ng ~ 515 in order to have a volume fraction 77 — 0.08 
corresponding to a porosity P = 92%. For the system 
simulated in the present work, system II, we fix the same 
77 and porosity increasing the box length to L = 20(75 . In 
this case the system consists of Ng = 1222 spheres. The 
average cluster size is Aug independent from the system 
size. The particles of the confined fiuid interact with a 
Lennard-Jones potential with a = ag aX variance with 
real systems where the aerogel spheres are much larger 
than the fluid particles. In the following Lennard-Jones 
units will be used. The potential is truncated at rc = 2.5. 
The interaction between the matrix spheres and the fluid 
particles is represented by a hard sphere potential with 
diameter a. The simulation of the phase diagram of the 
confined fluid is performed by Monte Carlo in the grand 
canonical ensemble (GCMC) with the algorithm intro- 
duced by Wilding^*-^^. By varying the chemical potential 
/i at constant temperature T we follow the behaviour of 
the density fluctuations and calculate the density distri- 
bution functions (DDF) P{p). In a two phase region the 
DDF develop a double peak structure. The coexistence 
point between the two phases is located with the Wilding 
criterionSS by tuning the chemical potential at constant 
temperature to obtain equal area under the two peaks. 
In the subcritical region where a large free energy barrier 
is present between the two phases the MESSLSSi22 has 
been implemented to enhance the sampling of the two 
coexisting phases. 

From the double peaked DDF we can extrapolate the 
biased sampling function at different chemical potentials 
and temperatures by means of the histogram reweighting 
techniques^. With the bias sampling function the MES 
can be carried on. In the final step the bias must be sub- 



tracted to recover the real distribution functionSa. To 
equilibrate system II lO'* Monte Carlo moves have been 
generated, while equilibrium properties have been calcu- 
lated with 10^ — 10^" Monte Carlo moves. The system 
has been investigated in the range of temperatures from 
T = 0.95 to T = 0.825. 



III. DENSITY DISTRIBUTION FUNCTIONS 
AND PHASE DIAGRAM 

In Fig.^the bimodal DDF of system II are reported in 
the thermodynamical range where the MES analysis has 
been performed. The low density peaks corresponding to 
the gas phase show a behaviour very similar to the corre- 
sponding quantities in system I therefore the gas branch 
of the coexistence curve appears as almost independent 
from the size of the system. On the liquid side coexisting 
with the gas we observe instead large differences in the 
structure of the peaks. In system I they were found to be 
broad and asymmetrio2&. Now in system II the form of 
the peaks is symmetric and sharper. To be more precise 
for temperatures T < 0.90 the density fluctuations are 
of the order of 0.07 in system I and 0.006 in system II 
in the same range of temperature variation AT ~ 0.075. 
By increasing the chemical potential at values above the 
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FIG. 1: Density distribution functions (DDF) of the fluid 
at the gas-liquid coexistence (solid line) for system II in the 
range of temperatures below the critical point from T — 0.825 
to T = 0.95 and liquid-liquid coexistence (long-dashed lines) 
in the range of temperatures from T — 0.825 to T = 0.85. 
DDF with closer peaks correspond to higher temperatures. 
In the inset P{p) of the liquid-gas coexistence in system I 
(solid line) compared with the liquid-gas (long-dashed line) 
and the liquid-liquid (dot-dashed line) distribution functions 
of system II at T = 0.85. 



liquid-gas coexistence we find in the adsorption isotherms 
the appearance of a second coexistence in system II. The 
corresponding DDF are also reported in Fig.^where they 
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show very well defined peaks. In the inset of Fig. ^ it is 
shown how the bimodal P{p) in the system I at T = 0.85 
is modified in system II with the emerging of new struc- 
tures. While the low density peak corresponding to the 
gas phase remains in the same position in system II, the 
liquid peak of system I splits in two in the larger sim- 
ulation box disclosing the presence of a possible LLC. 
The new bimodal structure is obtained by increasing of 
2% the value of the chemical potential from the value 
at the gas-liquid coexistence. The size of system I was 
not enough large to allow the resolution of the different 
contributions to the density fluctuations in the system in 
the region of the liquid peak. 

From the peak positions of the DDF with the Wilding 
criterion the coexistence curves of system II have been 
determined and reported in Fig. |21 where they are com- 
pared with the coexistence curves of system I and with 
the bulk. 
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FIG. 2: Phase diagram of system II (open circles) compared 
with the one of system I (stars) and a portion of the bulk 
liquid-gas coexistence curve (filled squares). In the inset 
liquid-gas coexistence curve of system II (open circles) with 
the diameter (open squares) and the best fit to the rectilinear 
diameter law (solid line), the filled diamond represents the 
critical point obtained by extrapolation of the fit. 



The confinement causes a substantial shrinkage of the 
gas-liquid coexistence curve and a shift toward lower tem- 
peratures and densities of the critical point. The confine- 
ment also does not induce large changes with respect to 
the gaseous phase of bulk fluid, in fact the gas branch 
of the gas-liquid coexistence curve in conflnement almost 
coincides with the analogous branch in the bulk. On the 
liquid side of system II a new coexistence curve is now ob- 
served between a low density liquid, which we shall call 
liquid I, and an high density liquid, liquid II. We note 
that the range of extension of the phase diagram of the 
conflned system is not quite affected by the system size. 
We checked that our results are independent from the 



matrix realization by performing simulations with other 
two different matrices generated with DLCA at the same 
porosity and L = 20. 

In system II, where the GLC appears as a well defined 
curve, it is possible to extrapolate the critical point po- 
sition with the use of the law of the rectilinear diameter 
and the scaling law of the densitj«2i. In the inset of Fig. El 
we show the coexistence diameter pd = {pg + pi)/2. The 
best fit to the rectilinear diameter law pd = Pc+A{Tc—T) 
combined with the scaling law pi— Pg = B{Tc — gives 
the following estimates of the liquid-gas critical param- 
eters: Tc — 0.96 and pc = 0.0791 to be compared with 
the values for the bulk = 1.1876 and pc = 0.3197iSa. 
At variance with experimental results^ we find a shift of 
the critical point to a lower density, this is due to the 
assumption of a pure repulsion between the fluid and the 
matri»i^. The exponent (3 compatible with these critical 
values is in the range 0.30 ~ 0.32 to be compared with 
the universal exponent (3 = 0.3258 found in the bulk^. 
This also confirms the experimental finding that the uni- 
versal behaviour of the gas-liquid transition in the bulk 
persists also in confinement^-. 



IV. LIQUID-LIQUID COEXISTENCE REGION 

In Fig. O we report the snapshots of two coexisting 
configurations of the fluid at T = 0.85 in the region of 
the second coexistence curve. From the snapshots we can 
observe that both liquid I and II are characterised by 
the presence of a large cluster and a number of smaller 
clusters with few particles. We calculated the gyration 





FIG. 3: Snapshots of the configurations of the confined fluid 
at T = 0.85 in the low density liquid phase (liquid I) on the 
left and in the high density liquid phase (liquid II) on the 
right. The gray and the black spheres represent the gel and 
the liquid particles respectively. 



radius Rg of the larger cluster in the two liquid phases 
in the range of temperature from T = 0.850 to T = 
0.825. Rg is deflned by Rg — VZ^i "^I J^^ where is 
the distance between the particle i and the center of mass 
of the cluster. In liquid I we have an increase of number of 
particles in the system upon increasing temperature and 
a 3% increase in the number of particles corresponds to a 
2.2% increase of the particles in the cluster and to a 7.3% 
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increase of the gyration radius. In liquid II we observe the 
opposite trend instead, an increase of number of particles 
is observed upon decreasing temperature. Besides a 5% 
increase in the number of particles corresponds to a 7% 
increase of the particles in the cluster while Rg grows 
only for an 1.7%. 

The radial distribution functions (RDF) of the fluid in 
the two liquid phases (Fig. ^ have been also calculated 
with a canonical MC starting from equilibrated configu- 
rations obtained from the GCMC. 
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FIG. 4: Radial distribution functions of the fluid in the two 
phases, liquid I on the left, liquid II on the right, for T=0.850 
(solid line), T=0.825 (long-dashed line). In the upper panels 
the RDF as obtained from simulation are shown together with 
the excluded volume correction (dashed line). In the lower 
panels the RDF are shown after the correction is applied. 



The usual algorithm for obtaining the RDF from sim- 
ulation consists in normalizing the number of atom pairs 
in a spherical shell with the uniform distribution of ideal 
particles in the same shell. Two of the grawi^) func- 
tions obtained with this procedure are reported as an 
example in the upper panels of Fig. ^ for the less dense 
liquid, liquid I, and for the denser liquid, liquid II. Ac- 
cording to previous literature^ when a fluid is confined 
the corrected RDF can be generated by appropriately 
normalizing grawif) with the uniform radial distribu- 
tion function guir) which depends on the geometry of 
the confining medium and a factor / = 1^/14// that ac- 
counts for the effective volume occupied by the particles: 
9{t) = graw{r)/{f ■ guir)). The function 5„(r) can be 
determined with a convolution of the structure factor of 
the aerogel system and the internal structure factor of 
the single hard sphere^. The resulting (?„(r) is reported 
in the upper panels of Fig.0] The factor / and 14// can 



be empirically determined by imposing that g{r) goes to 
1 for large r. For liquid I the factor / decreases and 14/ / 
increases with increasing temperature, in the two cases 
reported in Fig. 11/ = 2.03 for T = 0.850, TV = 1367 
and / = 2.18 for T = 0.825, iV = 1329. In the case of 
liquid II, in Fig. Hare reported T = 0.850, N = 1634 and 
T — 0.825, N — 1717, for which we get approximately the 
same value / = 1.81, independent from the temperature. 



V. CONCLUSIONS 

With the use of the DLCA algorithm to build a con- 
finig structure with fractal character we have found that 
the behaviour of the confined fluid confined in a fractal 
structure does not depend on the matrix realizations and 
finite size effects can be studied more systematically. By 
increasing the size of the simulation box the gas-liquid 
transition is well characterized, while a second coexis- 
tence curve appears on the high density side of the gas- 
liquid binodal curve. It is shown that in the low density 
liquid phase (liquid I) the size of the liquid domains in- 
crease with increasing temperature and equivalently the 
effective volume occupied by the particles 14// grows. In 
the high density liquid phase (liquid II) instead there is 
no increase of the size of the liquid domains and of the 
effective volume. In this regime we have an increase of 
the local density which is made possible by little rear- 
rangements of the particles. Even if drying effects due to 
the repulsive interaction between the fluid and the sub- 
strate cannot be excluded it seems that the high porous 
open structure of the aerogel network, built in our sim- 
ulation, allows the formation of an almost homogeneous 
liquid phase inside the free volume. 

The interpretation of fluid adsorption phenomena in 
aerogels is still controversial. In our simulation we stud- 
ied equilibrium phase transitions and we did not consider 
hysteretic behaviour. In this respect our treatment is 
analogous to the integral equation study of Krakoviack et 
al.l—, where a gas-liquid transition (but not a liquid- liquid 
one) is found for a fluid adsorbed in a high porosity aero- 
gel. Calculations performed with mean field density func- 
tional theory aimed to study nonequilibrium behavior 
show instead a different scenario. The disordered char- 
acter of the aerogel structure modifies the adsorption- 
desorption mechanism although at high temperatures 
there are indications of a gas-liquid transition^-^. We 
need future experimental and theoretical work to get a 
more complete understanding of this vaste and impor- 
tant phenomenology of fiuid adsoprtion phenomena. In 
this respect computer simulation with complete finite 
size scaling analysis seems to be necessary since we have 
shown here that finite size effects are quite relevant. 
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